{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys, inspect, time\n",
    "\n",
    "import numpy as np\n",
    "import torch \n",
    "import matplotlib.pyplot as plt\n",
    "torch.multiprocessing.set_sharing_strategy('file_system')\n",
    "\n",
    "import discrepancy, visualization\n",
    "from algorithms import ABC_algorithms, SMCABC, SMC2ABC, SNLABC, SNL2ABC\n",
    "from problems import problem_OP\n",
    "\n",
    "import utils_os, utils_math\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean-parma1 =  0.3720654037899112      mean-param2 =  0.21902819231692378\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAEGCAYAAAD2TVeiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1xUdf4/8Nd7hgHFQESxVFS0FEZoK6Uy081Lty2zi6137cK33dxES0vtO7mt2aTf7GJa/Xa9ZVmim7blJbdUMKVSRC3FxtrMLmol3vGCELx/fzDDgg4w4AxnLq/n4zEPmDNnZt6n6NXnnM/liKqCiIgqMxldABGRP2I4EhG5wXAkInKD4UhE5AbDkYjIjTCjC/BEs2bNNCEhwegyiCjIbN269ZCqxrl7LSDCMSEhAbm5uUaXQURBRkR+qOo1nlYTEbnBcCQicoPhSETkBsORiMgNhiMRkRs+C0cRmS8iB0Ukr8K2WBFZIyL/cf5s4u3vzcjIQEpKCsxmM1JSUpCRkeHtryCiEODLluMCALees20igHWq2gHAOudzr8nIyIDNZsOsWbNQWFiIWbNmwWazMSCJqNZ8Fo6qugHAkXM23wngTefvbwK4y5vfabfbYbfbsWLFCpw6dQq9evXCvHnzYLfbvfk1RBQC6vua48Wq+jMAOH82r2pHEfmTiOSKSG5+fr5HH+5wONC+fXu8/PLLWLx4MQCge/fucDgcXiidiEKJ33bIqOpsVU1V1dS4OLeze85jtVpx6tQpXH755Zg/fz4AIDs7G1ar1ZelElEQqu9w/FVEWgCA8+dBb364zWbD//zP/6BHjx7YsmUL5s+fj7S0NNhsNm9+DRGFgPqeW70cwH0Apjl/fuDNDx88eDAA4JlnngEAjB8/HrNmzSrfTkTkKfHVPWREJANATwDNAPwK4GkA7wP4J4A2AH4E8EdVPbfT5jypqala24Un+vfvj40bN2L//v2wWCy1rJ6IQoGIbFXVVHev+bK3erCqtlBVi6rGq+o8VT2sqn1UtYPzZ43BWFcPPvgg8vPzsWrVKl99BREFMb/tkLlQt9xyC1q0aFHeMUNEVBtBG45hYWEYMWIEPvzwQ/zyyy9Gl0NEASZowxEAHnjgAZSUlGDhwoVGl0JEASaowzExMRHdunXDG2+8AV91PBFRcArqcATKOmYcDgc2b95sdClEFECCPhwHDBiA8PBw3HbbbVyph4g8FhA32LoQK1euRHh4OAoLC3HkyBFs27YNaWlpAMDB4URUpaBvOdrtdkyZMgVnzpzBhx9+yJV6iMgjQR+ODocDI0eOxMUXX4wPPiibrciVeoioJkEfjlarFZ999hnuuOMOfPjhhzh79ixX6iGiGgV9ONpsNqSlpaFdu3YoKCjAK6+8wpV6iKhGQd8h4+p0efbZZwEA06ZNw2uvvcbOGCKqls9W5fGmuqzK484999yDnJwc/PTTTxARL1RGRIHMkFV5/NGdd96J/fv3Y+vWrUaXQkR+LqTCsW/fvjCZTOW91kREVQmpcGzatCl69OjBcCSiGoVUOAJlp9Y7d+7E3r17jS6FiPxYSIYjALYeiahaIReO7du3R0pKCsORiKoVcuEIlLUeN27ciCNHfHYLGyIKcCEbjiUlJbz5FhFVKSTDsUuXLoiJicGoUaO4xiMRuRX00wfdWbJkCUpKSlBcXIxjx44hNzeXazwSUSUh2XK02+2w2WwoLCzEJ598wjUeieg8IRmODocDo0aNQpMmTbBkyRIAXOORiCoLyXC0Wq3IycnBgAED8N577+HkyZNc45GIKgnJcHSt8dipUyecPn0aU6dO5RqPRFRJSHbIuDpdXNcYZ86cidmzZ7MzhojKhWTLESgLyLy8PNhsNpw+fRq9e/c2uiQi8iMhG44uQ4cORWlpKRYvXmx0KUTkR0I+HK1WK7p06YK3337b6FKIyI+EfDgCZa3H3Nxc7N692+hSiMhPMBwBDBo0CCaTCe+8847RpRCRn2A4AmjRogVuvPFGvP322wiEG44Rke8xHJ2GDRuG77//Hp999pnRpRCRH2A4Ot19990IDw/HHXfcwZV6iCg0B4G7s2LFClgsFhQXF+PEiRPIycnhSj1EIcyQlqOIPCYiu0QkT0QyRKSBEXVUZLfb8dRTT+HkyZP4+OOPuVIPUYiT+u6AEJFWALIBdFLVMyLyTwAfquqCqt6Tmpqqubm5Pq3LbDbj5MmT6NixI9q1a4cNGzaguLgYDRo0QElJiU+/m4iMISJbVTXV3WtGnVaHAWgoIsUAIgEcMKiOclarFZs2bcL48eMxevRobNiwASUlJVyphyhE1ftptaruB/ACgB8B/AzguKp+fO5+IvInEckVkdz8/Hyf1+VaqadDhw6Ii4vDuHHjuFIPUQir95ajiDQBcCeAdgCOAXhXRIapaqX5e6o6G8BsoOy02td1uTpdHn/8cRw6dAj5+fmYMmUKO2OIQpQRHTI3AtirqvmqWgzgPQDdDKjjPK6Veo4fP44mTZrA19c5ich/GRGOPwLoKiKRIiIA+gDwq/sTREVFYcyYMfjggw+wc+dOo8shIgMYcc1xM4ClALYB2OmsYXZ911GT9PR0REVF4bnnnjO6FCIygCHjHFX1aVVNUtUUVR2uqmeNqKM6sbGx+Mtf/oIlS5bgm2++MbocIqpnnD5YjbFjxyIsLAzXXnstpxQShRhOH6zGunXr0LBhQ5w4cQK7d+/Gvn37OKWQKESw5VgNu92OOXPmICwsDNOnT+eUQqIQUu/TB+uiPqYPumM2m1FYWIhx48bh9ddfx65du9C+fXtOKSQKEtVNH2TLsRpWqxXZ2dl46qmn0KBBA0yaNAnZ2dmcUkgUAhiO1XBNKdy1axceffRRvPvuuxg2bBinFBKFAHbIVMPV6ZKeng6HwwGz2YzY2Fh2xhCFALYca+CaUlhSUoIXXngBeXl5WLt2rdFlEZGPMRxrYeTIkWjTpg0mTpyI0tJSo8shIh9iONZCREQEnnnmGWzduhUJCQkcGE4UxHjNsZbCwsJgsVhQWlqKgoICbN68mQPDiYIQW461NHXqVDz99NPYv38/5syZw4HhREGKLcdacjgceOKJJ5CdnY2xY8ciNjYWgwYNgsPhV6uuEdEFYjjWktVqxaeffoply5ahX79+uO+++5CXl8eB4URBhqfVteQaGL5582a899576NKlC55//nk0b94cKSkp7KQhChJsOdbSuQPDk5KS0Lp1a2RlZeGxxx7D9u3bkZ2dzU4aogDHhSe8IDk5GVFRUdi8eTNmz56Nhx56CFlZWUhPT0deXp7R5RFRFfzxvtVBZffu3Thx4gTuvfdePPLII0hNTUX37t3ZSUMUwHjN0QusVitycnKwcOFCNGvWDEOGDMHatWvZSUMUwBiOXuDqpNm5cyfeeOMN7N69GwMHDuTqPUQBjKfVXnBuJ02zZs1w6NAhNGzY0ODKiKiu2HL0koqr9+zfvx9dunTBsGHDkJSUxOE9RAGI4egD4eHhGDJkCE6fPo3IyEicPn0as2bNgs1mY0ASBQiGo4/Mnz8f48aNw/bt2zFjxgzOwSYKMLzm6CMOhwPbtm3Djz/+CJvNhquvvho9evTg8B6iAMGWo4+45mDPnTsXiYmJGDhwIJYtW8bhPUQBguHoI67hPbm5ufjnP/+J06dP47777sMTTzxhdGlE5AGeVvvIucN74uPj8eOPPyIzMxMjRoyAiBhcIRFVhy1HH6o4vOeHH37A5MmT8dZbb6FVq1Yc3kPk59hyrEeXXnopIiMj8euvv2LNmjUQEa7eQ+Sn2HKsR1OnTsXixYvRoUMHDB48GB07duTwHiI/xSXL6pHZbEZhYSG+/fZbXHPNNUhJScHatWsRHR2NkpISo8sjCjnVLVnGlmM9slqtyM7OhtVqxfz587Fp0yYMHTqUw3uI/FCdwlFEwr1dSChwDe/JysrCXXfdhT/+8Y/44IMP8Pvf/97o0ojoHDWGo4isF5GECs+vAbDFhzUFrcGDB8NutyM9PR0NGjTAV199BavVirlz5+Kyyy5jDzaRH/Gkt3oqgH+LyEwArQD8AcADPq0qiA0ePLhSz/Tf//53jBo1CqdPn8bPP/+MXbt2sQebyA/U2HJU1Y8APAzgFQAPArhNVbf5urBQ8eqrr2LGjBk4dOgQHnroIS5QQeQnamw5isgkAAMA/B7A7wCsF5Fxqrqqrl8qIjEA5gJIAaAAHlTVz+v6eYHM4XDgz3/+M06dOoWJEydi48aNvP8MkR/wpEOmGYBrVPVzVf0HgFsAPHqB3/sKgH+rahKAKwCEbBK4erDT09PRvHlzPPPMM+U92kRkHE9Oq8eo6pkKz39Q1Zvq+oUiEo2yVug85+cVqeqxun5eoHP1YG/evBljx47F2rVrMWzYMN5/hshgRkwfbA8gH8AbInIFgK0AxqjqKQNqMdy5C1SYzWY0a9aMnTFEBjNiEHgYgM4A/p+qXgXgFICJ5+4kIn8SkVwRyc3Pz6/vGutVxQUqnn32WezYsQM5OTlGl0UU0owIx30A9qnqZufzpSgLy0pUdbaqpqpqalxcXL0WaKRHHnkEsbGxmDJlitGlEIU0TwaBNxCRR0TkdRGZ73rU9QtV9RcAP4lIonNTHwBf1fXzgk1UVBQee+wxrFy5Etu2ccQUkVE8aTkuBHAJynqpPwEQD6DgAr83HcA7IrIDwJUAnrvAzwsq6enpaNiwIXr27MlZM0QG8aRD5jJV/aOI3Kmqb4rIIgAfXciXquoXANyuhEHAhx9+iIiICBw7dgxbtmxBQUEBZ80Q1TNPWo7Fzp/HRCQFQGMACT6riGC32/Hmm28iKioK06ZN46wZIgN4Eo6zRaQJgKcALEfZ9cH/82lVIc7hcOAPf/gDRo8ejWXLlmHnzp2cNUNUzzwJx3WqelRVN6hqe1VtDuBjXxcWylyzZsaOHYuoqChMnjyZs2aI6pkn4bjMzbal3i6E/ss1a+bLL79Eeno6li1bhuHDh3PWDFE9qrJDRkSSACQDaCwi91R4KRpAA18XFsrOnTVjMpnQokULdsYQ1aPqWo6JAPoCiAFwR4VHZwAP+b600FZx1sxf//pX5ObmctwjUT2q8QZbInKd0cuJBcsNturq+PHjSEhIQI8ePbB8+XKjyyEKGhd6g63t3pwhQ7XXuHFjjBs3DitWrEAo/0+CqD4ZNUOGamn06NFo1KgRevXqxVkzRPXAkBkyVHurVq1CeHg4jh49iuzsbBQVFXHWDJEPcYZMgLDb7Vi4cCGaNm2KZ599lrNmiHysNjNkJuG/M2Se92lVdB6Hw4Gbb74Z48aNw7///W84HA7OmiHyIU9ukzDXOUPmE9cMGVX9e30UR//lmjXzwAMPwGQy4e233+asGSIfqm4Q+Njq3qiqL3m/HKqKa9bMvHnzcOONN2Lu3LnIyMjgaTWRj1TXIRPl/JkI4GqUnVIDZQPBN/iyKDpfxVkzX331FVQVkyZNYmcMkY9UeVqtqpNVdTLKbs3aWVXHqeo4AF1QNpyH6plr1kxBQQEaNWqEAwcOGF0SUdDypEOmDYCiCs+LwN5qQzVq1Aj9+/fHu+++i8LCQqPLIQpKng4CzxGRv4nI0wA2A3jTt2VRTYYPH44TJ05gxYoVRpdCFJQ86a22A3gAwFEAxwA8oKpTfV0YVa9Xr15o2bIlFi5caHQpREHJkxkyUNVtALgkjB8xm80YOnQoXn75ZeTn5yOUbl9LVB+MuG81ecnw4cPx22+/YcmSJUaXQhR0GI4B7PLLL8cVV1zBU2siH2A4BriUlBTk5ORwpR4iL/PomiP5p4yMDGzcuBEiggkTJuCmm27iSj1EXsKWYwCz2+1YsGABbrzxRmRkZOCGG27gSj1EXsJwDGCulXmGDRuG77//Hrm5uVyph8hLeFodwFwr9dx+++0wmUxYtWoVTp06xZV6iLyALccA5lqpZ8eOHbj22muRkZGBtLQ03t+ayAvYcgxg7lbqee2119gZQ+QFbDkGONdKPTt27AAAhIeHG1wRUXBgOAaJ5ORktGnTBitXrjS6FKKgwHAMEiKCvn37Ys2aNVzGjMgLGI5BpG/fvjh9+jQ++eQTo0shCngMxyDSq1cvREZG8tSayAsYjkGkQYMG6NOnD1auXAlVNbocooDGcAwyffv2xffff4+vvvrK6FKIAhrDMcjcfvvtAIBVq1YZXAlRYGM4BplWrVrhqquucnvdMSMjAykpKVzejMgDhoWjiJhFZLuIsPfAy26//XZ8+umnOHLkSPm2jIwM2Gw2zJo1C4WFhZg1axZsNpvbgGSIEhk7fXAMAAeAaANrCErh4eEoLS1Fs2bN0KlTJ9hsNtjtdkyZMgU7d+5EVFQUevbsiXnz5mHEiBGw2+1wOBywWq3o1asXVq1ahXnz5qF79+7Izs7mGpEUmlS13h8A4gGsA9AbwMqa9u/SpYuSZxYtWqTt2rXTxo0b66BBg3TdunXaokULBVDp0bFjR7377rsVgGZmZmpRUZGuW7dOLRaLPvjgg7p8+XI9fPiwqqpmZmZqcnKywUdG5H0AcrWqnKrqBV8+ACwF0AVAz6rCEcCfAOQCyG3Tpo3v/ukEmeTkZM3MzNT7779fo6OjtUuXLuWBOGTIEP3qq690zpw52rNnz/LtV1xxhSYnJ2ujRo0qBWh0dLQ+/fTTevDgQRURTU5OVpPJpMnJybpo0SJdtGjReduIAolfhSOAvgBed/5eZThWfLDl6DmTyaRFRUW6fPlyBaBJSUn62muvqYhou3btyluJmZmZCkAHDRqkPXv21DvvvFMfffRRveSSS1RENDMzU/v3768ANCIiQk0mk65cubL8vXFxcRoXF1fp89q1a6ejRo1iYFLA8LdwnApgH4DvAfwC4DSAt6t7D8PRc66WY2lpqX799ddaUlJSflp8bksvPj5eMzMzK73fZrOpxWIpD73Zs2eriCgAbdq0qd588816zz33aKNGjTQqKkonTJigb775pu7du/e897oCkwFJ/sqvwrHSl7Pl6HWua46eBFRV+57b+hMR/fzzz3XAgAF67bXXanJycvmpt8ViKf89LCxMAejcuXP1yJEjqsrrleTfGI4hpjbXAj3Z19UarSghIUETEhK0pKREv/zyS501a5YCUJPJpAD08ssv15MnT2pRUZGaTCavHyORN/htOHr6YDgay10L0901R4vFov/7v/+r77//voqIDhkyRNetW8eWI/kthiNdMHctzHO3jRo1qjxEp0yZogA0NjbW69ccPamF1znJEwxHqjcVQyoqKkpFRNevX+/Vz/ekFcuOIPJEdeHIudXkVa572pSUlGDfvn3o2LEj+vXrh8TExBqnI3oybdFut2PevHm4/vrrkZmZiYKCAogIAODgwYPYsmULevXqhXnz5sFut/v0WCnIVZWa/vRgyzFwTZ8+XUVEk5KStKCgoMpWXcUW4dmzZ6vczzWOc+DAgefN+nE9Vq1axY4g8gh4Wk1GSU5O1r/97W8KQAcOHKhnzpxxO7wnOTlZ33//fb3lllu0c+fOevLkySr3e+655xSAPv7445qbm6stW7bUli1bal5enqakpGirVq10+fLl7AiiGjEcyTCult60adMUgKampuqePXvOa9WZTCa97LLLNDw8XEVE77vvPretv/nz56vZbNa2bduWB2jFa46ff/65mkwmveiii3jNkWpUXTgauSoPhQCr1Yrs7GxMmDABSUlJGD58ODp37ozY2FikpKTA4XCgdevWUFXk5+cjMzMTa9asweTJk3HxxRfDarVW+jzX9UyTyYTo6GhYrVa88sorAID09HQ4HA7Exsbi0KFDiI2NNeKQKVhUlZr+9GDLMXCd27u8YMECNZvNCkBHjRqlM2fOVJPJpGazWZs2baqZmZl65swZvfLKK1VE9Pnnny//rNzcXDWZTPrwww9X+51nzpxRq9Wq8fHxeuzYMV8fIgUw8LSajHTuGMSWLVtqt27dyjtQ+vbtqytXrtT4+Pjy/RITEzU6Olo7deqkJ0+e1OLiYu3cubNecsklevTo0Rq/c9OmTSoi2qRJE459pCpVF448rSafGzx4cKWFcs1mM7777ju8+uqrOHPmDJ588kmUlpbiwIED+Omnn8r3W7t2LW666SbEx8fj+PHjUFWMGTMGMTExNX7nd999h+joaBw9ehSrV69GREQEF+2l2qkqNf3pwZZjcHE3V9tdz/SiRYs0JiZGAajZbNauXbt6PLg7OTlZV69erYmJidqmTRv97rvvuAgGnQc8rSZ/4unKQcnJybpmzRq94YYb9KKLLtIffvjB44Bz9ZJv3rxZGzdurNHR0bpw4UKOfaRKGI7kdzyZC+0KuMLCQt2/f7+qqseDuyu2Tvfu3avXXXedAtDIyEi1Wq28DkmqynCkAOXp6bc757ZO16xZo5GRkQpA4+PjddOmTZyDTQxHCky1Wbi3qvefu/L5Sy+9pK1atdLw8HBdv349r0OGuOrCUcpe92+pqamam5trdBlkgIyMjEq3jrXZbHXubTabzSgsLMSJEydw/fXX4/jx48jJyUFCQgJKSkq8XDkFAhHZqqqpbl9jOFKoSElJwaxZs9CrVy/k5eXhmmuuQWJiIoqLi5GXl2d0eWSA6sKRS5ZRyLDZbEhLS0NWVhYSExMxevRofPHFF+jQocMFfa4nS61RAKrqfNufHrzmSN5y7nVI1/27x48fX6fVxS/0uigZC+yQIXLv9OnT2qZNGzWZTJqRkVHr+3LHx8fr888/r/369dOlS5eqKu+4GEiqC0dOH6SQ1rBhQ0RERCAiIgIvvfQSDh8+jKNHj+Ls2bMoKSnBP/7xD2zatAl9+vTBoEGD8MILL+Cjjz5Ct27dMH36dEyaNAnjx4+H2WzGRx99hNatW6N79+5wOBxGHxpdIIYjhbw9e/Zg0aJFGDp0KLZs2VLptZycHCxZsgQAYDKZUFpais8++wyPPvooduzYgbCwMDRu3Bg7duzA9ddfj7vvvhszZsw4b6k1CjzskKGQZ7Va0bx5c/z66684ePAgzp49i4SEBCQkJOC7777Dr7/+ioyMDJSWliIsLAxPPfUUioqKsHDhQsydOxeHDx/G119/jWXLluHIkSMYPnw4nnjiCaMPiy4Qw5FCnqsXe8eOHYiJicGnn36KU6dO4dSpU8jKykKTJk1w8cUXw2KxYMKECfjll1+Ql5eHYcOGoU2bNoiPj0d6ejquvvpqxMXF4ezZs5g7dy57sANdVRcj/enBDhnytdrel7u6nul7771XAeif//xnr/Vg877cvgH2VhN5hychlZycrDfccIOKiH700UeqemE92J7emZFqr7pw5AwZIi8zm804fPgwrrvuOogIdu7cidLSUjRo0KBO0xRTUlIwc+ZMZGRk4IsvvsD69euRk5OD9PR0zuy5QJwhQ1SPrFYrtm/fjsmTJ8PhcGDZsmXIzs6ucw+2w+HAl19+iblz5yI3Nxfjxo3jcKF6wKE8RF7m6uCZPXs2EhMTMXHiRACA3W6v0+e1bdsWjz/+OPr164cOHTrgxRdfRHx8PIcL+VpV59v+9OA1Rwo0rmuTIqIA9LHHHqvT5xw4cEAbN26sYWFhunz5ci0oKNBLL71UTSaTvv7667WqhZ055wM7ZIiMUVxcrB06dNArr7xSS0tLPXpPxTCLjIxUi8Wi06ZNK9922WWXqcVi0VtvvbXGz+Tc7+oxHIkMtGDBAgWgy5cvr3HfimE2evRoBaBxcXHnhdmrr76qALRFixY19pxnZmbqwYMHdfHixVpcXMy53xUwHIkMVFRUpO3bt9fU1NTzWnruVitfvHixPvfccwpAR40a5TbM3nnnHW3YsKFaLBbdvn17lS1Ck8mkO3fu1Hbt2ikAvfbaa3XXrl280ZgTw5HIYHPnzlUA2rZt2/MGlX/88ce6bt06HTRokAIof/Tu3VvPnj3r9qZiycnJunTpUo2Li9MOHTrohg0b3IZou3btNCoqSuPi4nTatGkaExOjERERGhMTU6cl2oINw5HIYG+99ZaGhYVpUlKSnj17VteuXathYWF61VVXadOmTRWAhoWFaYMGDbRJkyaal5dX3sp0F3quOzOuX79e4+PjFYDefffdKiLlAde6dWs1mUwaFhamCxcu1KKiIl2yZIlaLBYFoF27dtWffvqp2iXagj0gGY5EBktOTtbHHntMAejQoUPLAw2ADho0SJcuXarHjx/Xt956SwF4dE9v150ZT506pc8884yGhYUpAB0wYIDa7XYFoBEREZqWlnbeqfsjjzyiERERGhkZqfHx8Wo2m9VsNmuHDh20S5cu5aEZ7NcmGY5EBjOZTFpQUKAJCQlqsVj0jjvu0FatWqmIVNovMzNT4+Pj67QCucVi0d/97nflw4f69++vq1evrrLVmZeXpyNHjtT777+//D0DBgxQAPriiy96fI/wQOZX4QigNYAsAA4AuwCMqek9DEcKdK6W3uHDh/XIkSOqqmqz2dRisXjt1rMiokVFRbpt2zadM2eOlpSUVHm98tz7gSckJGhCQoKqql5++eXaq1cvn7UcvX1t80I+z9/CsQWAzs7fowB8A6BTde9hOFKgq2q84bm3XbiQoHAXeu4Czl0tFa85jh8/Xk0mk7Zt29br1xy9Pe7yQj/Pr8LxvAKADwDcVN0+DEcKBr7uDa5NUNS0RBsATU9P92p9qv8N8KlTp+qCBQu0tLS0yhZqTf+8Dh8+rO3atdOxY8fqmDFj9OGHH1bV2q2A5LfhCCABwI8Aot289icAuQBy27Rp49GBEoU6bwTwb7/9ps2aNdNhw4Z5vT6TyaTvv/9+eWfUwIEDNT8//7xT/4pBf+DAAX355Ze1adOmevPNN2vv3r31kksuqTTsKTIyUnv06KGlpaW1ulbql+EI4CIAWwHcU9O+bDkS1a8RI0ZobGysFhcXe/VzO3XqpC1atNCkpCR99tln1Ww2a4sWLbR58+aVQr1ly5Y6cuRITUpKqhSCJpNJr7nmGn3ggQd0+vTp5af+JSUl5d8R0C1HABYAHwEY68n+DEei+vXuu+8qAN24caPH7/Gk1epaJf2FF17QoqIinTlzZnlP+ciRIzUzM1P79OlTHoZdu9OGyPEAAAgoSURBVHbV6dOn6+rVq3XPnj3VtjAD/pojAAHwFoAZnr6H4UhUv44fP65hYWE6YcIEj/b3JKT27NmjDRo00K5du57XSuzWrVt5IEZHR2tUVJS2bNmy0nfU9dpkdfwtHLs7/yHsAPCF83Fbde9hOBLVvz59+minTp082jc5OVlXrVqlaWlp+s4777jtaOnbt682atRI9+3bV+m9JpNJz549qxkZGfrGG2/oyZMnPR4Mf6H8Khzr8mA4EtW/GTNmKADds2dPjfuaTCbt169feeuvf//+un///krTGQHokCFDzntvVUOQPBkMf6EYjkRUa99++60C0JkzZ9a4b7NmzcqvJU6bNk0tFotGRkaq2WzW1atXa9u2bbVt27ZuW39GrjnJcCSiOklKStKbbrqp2n1c61VGRUXpunXrtKioSOfMmVPe0ZKSkqIANCsryyfXDS8Ew5GI6uSJJ55Qi8WiJ06cqLT93AHjycnJ+tZbb503nXHChAlqMpl06NChqqp+N1+7unDk3QeJqEp9+/ZFcXExPv744/JtGRkZsNlsePLJJ9GkSRPEx8fj5MmTCAsLQ15eHkpKSpCXl4dOnTrhlltuwd69ezF//nwAuKC7MNa7qlLTnx5sORIZo7i4WBs1aqQxMTHlLcJWrVrppEmTtGPHjtqkSRP95ptvPJ7D7W9rRIKn1URUF4sWLdJGjRpp48aNNSsrS/v27VveI92qVStdv369qlZ9uuzvq4tXF45S9rp/S01N1dzcXKPLIAo5KSkpuOuuu8rvuR0ZGQkRwUUXXYT9+/fDbDYDALKyspCeno68vDwjy601EdmqqqnuXgur72KIKHA4HA58+umnKCgoQOfOndG/f3/861//wogRI7BhwwZ0794d2dnZSEtLKw/QYMFwJKIqWa1WbNu2Da+88kr5tvj4eMTHxyM9PR0OhwNWqxV2ux2DBw82sFLvYzgSUZVsNhvS0tIwb968Sq3E559/PujC8FwMRyKqkisAg72V6A47ZIgoZFXXIcNB4EREbjAciYjcYDgSEbnBcCQicoPhSETkRkD0VotIPoAfavm2ZgAO+aCc+hYsxwHwWPxRsBwHULdjaauqce5eCIhwrAsRya2qiz6QBMtxADwWfxQsxwF4/1h4Wk1E5AbDkYjIjWAOx9lGF+AlwXIcAI/FHwXLcQBePpagveZIRHQhgrnlSERUZwxHIiI3gi4cReRWEflaRL4VkYlG11MbIjJfRA6KSF6FbbEiskZE/uP82cTIGj0hIq1FJEtEHCKyS0TGOLcH4rE0EJEcEfnSeSyTndvbichm57EsEZFwo2v1lIiYRWS7iKx0Pg/IYxGR70Vkp4h8ISK5zm1e+xsLqnAUETOA1wD8AUAnAINFpJOxVdXKAgC3nrNtIoB1qtoBwDrnc3/3G4BxqmoF0BXAI85/D4F4LGcB9FbVKwBcCeBWEekK4P8AvOw8lqMA0gyssbbGAHBUeB7Ix9JLVa+sML7Ra39jQRWOAK4B8K2qfqeqRQAWA7jT4Jo8pqobABw5Z/OdAN50/v4mgLvqtag6UNWfVXWb8/cClP2H2AqBeSyqqiedTy3OhwLoDWCpc3tAHAsAiEg8gNsBzHU+FwTosVTBa39jwRaOrQD8VOH5Pue2QHaxqv4MlIUOgOYG11MrIpIA4CoAmxGgx+I8Df0CwEEAawDsAXBMVX9z7hJIf2czAIwHUOp83hSBeywK4GMR2Soif3Ju89rfWLDdJkHcbONYJYOIyEUAlgF4VFVPlDVSAo+qlgC4UkRiAPwLgNXdbvVbVe2JSF8AB1V1q4j0dG12s6vfH4vT9ap6QESaA1gjIru9+eHB1nLcB6B1hefxAA4YVIu3/CoiLQDA+fOgwfV4REQsKAvGd1T1PefmgDwWF1U9BmA9yq6jxoiIq3ERKH9n1wPoJyLfo+ySU2+UtSQD8VigqgecPw+i7H9a18CLf2PBFo5bAHRw9r6FAxgEYLnBNV2o5QDuc/5+H4APDKzFI87rWPMAOFT1pQovBeKxxDlbjBCRhgBuRNk11CwA9zp3C4hjUdUnVTVeVRNQ9t9GpqoORQAei4g0EpEo1+8AbgaQB2/+jalqUD0A3AbgG5RdF7IZXU8ta88A8DOAYpS1gtNQdk1oHYD/OH/GGl2nB8fRHWWnZjsAfOF83Bagx/I7ANudx5IH4K/O7e0B5AD4FsC7ACKMrrWWx9UTwMpAPRZnzV86H7tc/61782+M0weJiNwIttNqIiKvYDgSEbnBcCQicoPhSETkBsORiMgNhiP5JRGJEZG/VHjeUkSWVveeOn5PTxHp5u3PpcDHcCR/FQOgPBxV9YCq3lvN/nXVEwDDkc7DcY7kl0TEtaLS1yhb7OE1lA1aThGR+1G22ooZQAqAFwGEAxiOsiXGblPVIyJyqfN9cQBOA3hIVXdX+I4EAJsAlADIB5Cuqhvr4/jI/wXbwhMUPCYCSFHVK4HyIKsoBWWr/TRA2cyOCap6lYi8DGAEyuYMzwbwsKr+R0SuBfA6yuYTAwBU9XsR+TuAk6r6go+PhwIMw5ECVZaWrRVZICLHAaxwbt8J4HfOFYG6AXi3wmpAEfVfJgUqhiMFqrMVfi+t8LwUZX/XJpStU3hlfRdGwYEdMuSvCgBE1fXNqnoCwF4R+SNQtlKQiFzh7e+h4MVwJL+kqocBfCoieSIyvY4fMxRAmoi4Vm5xd8uMFQDudt6kqUcdv4eCEHuriYjcYMuRiMgNhiMRkRsMRyIiNxiORERuMByJiNxgOBIRucFwJCJy4/8DkpZkZ8CyevMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAFBCAYAAAAsfIegAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3hU1daH351eCAmht9CrWFGKAmIXUdQrtmu7luv181qvBdu1947lYkWxV0AEEYHQpEoTBAIJIZWSXifJtPX9sSd9JpmESYHs93nOc2Zy9jmzz2TmN2uvvfZaSkQwGAwGg+/wa+kOGAwGw9GGEVaDwWDwMUZYDQaDwccYYTUYDAYfY4TVYDAYfIwRVoPBYPAxAS3dgaamU6dO0rdv35buhsFgOMrYtGlTloh0dnfsqBfWvn37snHjxpbuhsFgOMpQSiV7OmZcAQaDweBjjLAaDAaDjzHCajAYDD7GCKvBYDD4GCOsBoPB4GOMsBoMBoOPMcJqMBgMPsYIq8FgMPgYI6wGg8HgY4ywGgwGg49pNcKqlOqtlFqmlNqllNqhlLrbTRullHpLKZWglNqmlDqpJfpqMBgMddGacgXYgftEZLNSKgLYpJRaLCI7q7SZBAxybaOBGa69wWAwtBpajcUqIgdEZLPrcSGwC+hZo9nFwGeiWQdEKaW6N3NXDQaDoU5ajbBWRSnVFzgRWF/jUE8gtcrzNGqLr8FgMLQorU5YlVLtgB+Be0SkoOZhN6fUqt+tlLpVKbVRKbUxMzOzKbppMBgMHmlVwqqUCkSL6pciMttNkzSgd5XnvYD9NRuJyAcicrKInNy5s9s8tAaDwdBktBphVUop4GNgl4i87qHZPOB6V3TAGCBfRA40WycNBoPBC1pTVMBpwHXAdqXUVtffHgFiAETkPeAX4AIgAbAAN7ZAPw0Gg6FOWo2wisjvuPehVm0jwL+bp0cGg8HQOFqNK8BgMBiOFoywGgwGg48xwmowGAw+xgirwWAw+BgjrAaDweBjjLAaDAaDjzHCajAYDD7GCKvBYDD4GCOsBoPB4GOMsBoMBoOPMcJqMBgMPsYIq8FgMPgYI6wGg8HgY4ywGgwGg48xwmowGAw+xgirwWAw+BgjrAaDweBjjLAaDAaDjzHCajAYDD7GCKvBYDD4GCOsBoPB4GOMsBoMBoOPMcJqMBgMPsYIq8FgMPgYI6wGg8HgY1qNsCqlZiqlMpRSf3k4PlEpla+U2uraHm/uPhoMBoM3BLR0B6rwKfAO8FkdbVaJyIXN0x2DwWBoHK3GYhWRlUBOS/fDYDAYDpdWI6xeMlYp9adSaqFS6piW7ozBYDC4ozW5AupjM9BHRIqUUhcAc4FB7hoqpW4FbgWIiYlpvh4aDAYDR5DFKiIFIlLkevwLEKiU6uSh7QcicrKInNy5c+dm7afBYDAcMcKqlOqmlFKux6PQfc9u2V4ZDAZDbVqNK0Ap9TUwEeiklEoDngACAUTkPWAq8H9KKTtQAlwlItJC3TUYDAaPtBphFZGr6zn+Djocy2AwGFo1R4wrwGAwGI4UjLAaDAaDjzHCajAYDD7GCKvBYDD4GCOsBoPB4GOMsBoMBoOPMcJqMBgMPsYIq8FgMPgYI6wGg8HgY4ywGgwGg48xwmowGAw+xgirwWAw+BgjrAaDweBjjLAaDAaDjzHCajAYDD7GCKvBYDD4GCOsBoPB4GOMsBoMBoOPMcJqMBgMPsYIq8FgMPgYI6wGg8HgY4ywGgwGg48xwmowGAw+xgirwWAw+JiAlu6AwWA4fOx2KCwEiwVKS6GsrPZWWgo2GzgcenM6Kx+XPxcBPz9QSu9rPg4IgKAgCA7WW83H4eEQEaG3gDasLq3q1pVSM4ELgQwRGeHmuAKmAxcAFuAfIrK5eXtpMDQNxcVw4ABkZ3ve8vK0SNYkIADatYOwMAgJqRS7qltICAQGgr9/5ebnp/fBwXoPWlzLRbZcfO32SgEuF2qrtfbj4mIt8IWF+pyqKAUdOkCXLtC1q96Xb+XPo6N1uyOdViWswKfAO8BnHo5PAga5ttHADNfeYGjVOJ2QlAR79kBamt7S0/XeYtFtwsKge3fo1Ak6dtRbnz5w0kmVzzt00CJ4JOJ06h+GQ4cgI6Ny++svvT90SP94iGihHTgQBgzQ+4EDoWdP/UNwJNCqhFVEViql+tbR5GLgMxERYJ1SKkop1V1EDjRLBw2GerBaIT4edu6EXbv0lpamBaFPHxg8GGJiYMIE6NVLi0V4eEv3unnw89MWaXQ0DBvmuZ0IZGZCQoLeli2DDz+E1FRtVZ90Eowdq7devZqv/w2hVQmrF/QEUqs8T3P9zQiroUVITYU1a/S2dav+4g8erIVjzBi48Ub95W+q4W2JFQ4VwKF8yCiE7CIoKoOi0ur74jK9L7GCwzXEd9cnEf13pSA4oMoWqPftgqFDOESHV9mHVT7v2A78D9OqVKrSRXDqqdWPlZXB5s2wdi189522+nv31iJ76qladFuDb7cVdKFBuPt4Sq1GSt0K3AoQExPT1H0ytAHy8/Uwfs8ePaT/809ticbE6C/0ddfBq69qH+bhUFymxbF8yymG7GL9OKMA0nL138T1qQ8Ngi4R0DUSurbXAte1PQzoDO1CtBBW7IMhJFALX7l4esLpBKsDymxQZndtNi3OOcWQW6z323Mrn5f30yn6NXpEQZ+O0L9z5RYTfXjD+eDgSmu1nNRULbRffQX/+Y8eBTz9NAwZ0vjXOVyUSC1dalFcroD5Hiav3geWi8jXrue7gYl1uQJOPvlk2bhxYxP11nA04XTqYXxcnBbQ3bshMVHPpEdG6i/qoEHQv7+2SHv39v7aNjuk5EBSFuzPgwP5cMC1P1SgXxsgPFhbfR3DIbp877IEO0dATEdtIbZ2X6Pdoe8tOQsSM/WWkAH7sqB7JEwcChOHwPAevr+XnTvhppvgm2+gb1/fXrsqSqlNInKyu2NHmsU6D7hDKfUNetIq3/hXDY0lNxfWr4d162DDBj2TPWgQDB8OQ4fClClaRL21QstsEHdAC0hiJux1CUqpDQL9tbXWr7O25I7pAWcP1yLTpb374bMI5DvhkB0y7LDXAeutUFAKBU7X5oBC12Ob6OFbVUO03Gzy1hNR3t4PaO8HUf4Q6Q9RfnrfPQD6B8KAIOjk79nqDfCH3tF6Gze4+rH0XFixG95aAjv2VxfaY3oevttk+HCYOROuvhpmz9YTgs1Nq7JYlVJfAxOBTsAh4AkgEEBE3nOFW70DnI8Ot7pRROo0R43FaignKwvmzYPff9fWaFQUjB6th5WnnKKfe0thCfyZCltS9BZ/SIvnsO4wqCsM6KKHvv06QZiHWfw8ByRaIdEG+1z7FBsUOSvbRPlB1wDoEqD3nfwh0g/a+0OEnxa/9n4Q4Q9BPvTjOkQLdr4D8qrs99t0PxOtkOHQbTv5a7HtHwQjQ2BkKAQ0oC/lQrt4h7bqH5kMZw47fIHdtAnuvlv/z6OjD+9a7qjLYm1VwtoUGGFt2+TlwZw58MMPOq7y0kth4kQ9wdSQIWh6Liz6C5bshORsiAiB43vDSX3gxBgY2MX99UQgzQ5bSmFrKWwvrRSkKD/oFwT9XKLULxD6BGqRrIkDIY9S8ijFgg0Ldte+6manFDuC4EQQ13mVzwWFIoQAgvEnhIBqj4MJoD1BdCSUToQRRTCqHltXBHIcWmwTrLChBDaWQLQ/nBEOZ4XDMcHg56VIpuXA8wtgbwY8dhGMH1z/OXWxciU88YQW14iIw7tWTYywGmFtUxQU6C/S99/rGNFLLoHLLoNu3by/RpkNVifAr9th7V49XD9/BJxzjJ6QcWdNOQV2W7WIbimFbaVgcULvQDghBE4MgWNDoKubIXQRVlIpYD9FZFNCDiVkU0IhVhTghyKKEKIIIZxAwggkjADXvnILxh8/VMVWfm755kAow04pDte+8nEZDgooIwsLmZRQQBmC4I8f0YTQiTB6EcEgoulJBH51iG6GHZYVQ2wx/FWm34Mzw+HiCG1510dyFjw7X/+gPT4Fxgzw/n9Xk4ULYfp0mDtXL5LwFUZYjbAe9VitsHgxfPqpDjKfMgUuv1zPEHtLUpYW0oXbIdcC4wZpMR07AALdiIFNtAW6ygIrLFpMhgbr4XC5iIZXsWIFIZsS0ikklULSKCCdQmw4CSeQ3rSnB+3oSBgdCaUjobQjsJbV6MSGlWKXnVqMjaKKxw6srlaKcs+qqvLcn0ACCSeUjkTQg0DqD6K14ySHErIoIYV84sklnUIC8aM/HRhCNGPoWafQJlthaTF8XwBhfnB3NEzwIn53bwY8PU//P164TPtgG8P338OsWfoH11eTZUZYjbAelRQWamtk7lwdcnPGGXD99XqVjjdkFsCyOD28356uLdHzjoFJx0G3yNrt022w3KKHu9tK9d9GBMO4MJgQBt2rTHKV4SCFfPaRRxL5pFKADScdCaUnEfSmPb2IoCcRBKHH/oKDEnKwkOWyVyu3UvIQtPNVi2MYgbQjsMJ+1Xt/gmr0WlxOAf09d2DDSiElZFPEfmxYUPgTThci6Ek3TiKK/l69f2XY2UsesSTRnw5cgHdm5Z4yuHE/rO7nVXMAHvoeekXDHWd5f05V5s6FF16AVat0TgNfcDRFBRjaOJmZ8NNPerNYYNIkHbPojZjaHRC7S/tKN+yDqDA4Y6j+so5ws1wy266FdGkx/FmqZ8QnhsP1kXBs18rJIgdOksjnZ7KII5tsSgjCnxja048ozqIvMbQnEH8EoZRcithPIbuII51C9mOlEIWfy07tRBgdCacLnRhOKNGEEIXCjfPVBzixU0wGRaSzix/ww58RXEM4dftOgglgOJ0YRDSPspwx9CCa0Hpfzyb6B8lb9mbAukRYepn351Tl44/h55/1Ci5fiWp9GIvVcESQmAivvKLjS6+6Ci66yHuf6fY0mLUafo/Xs82Tj4NR/WoP722ifYKLi2BDqZ5tnximJ2COD6mcgCnDzm5y2OkS0lLs9CWS4XRiGJ3oRCgKheCgkHRy2UsuCeSTjBM7IXQggp60owcR9KAdPQimfWVHxArOTHBm6L0jw/U4AyQXxA7YK/fVHgeBfzfwc23+3cCvu2vfFVT9ipbNbnbwJe3pzTCurN43D+wkiwUk8ABj6m37Zb6OMrjdi5l6EbjgDXjpcjiuAXHD5ee+9JJeVvzRR4e/eKMmxmI1HLFs366/HPn58MADMH68d2E4mQXw9Xr4cRP07QQ3nAYvX17bKhWBP0rh8zy9PzscLm0Pz3etHr6UTiFrSWcLB1EoBhPNMXTiQgbSzjX8tpBFJuvZSgL5JCE4iaAnHRhIDBOJpE/lUN2ZDfZ4sP8Bjj36sXO/69WCwK+La+sM/l0gYAD4jQUVBSoICAAVUH1PAFAGjoPgdG32BHD+rh87DoIUQ9BECL0cAk5w+2Z2ZAjjeYqDbGQ1z9GdkQziYgLwLMrD6cRSkthOBsfSpc7/zZ+lcImXM/Sf/A4n9224qDqd+vMiAp980vwLKozFamiVrFkDL7+sZ3GnTYMTT6z/HIcTFvwJn64GixWuHg2XjdRLOmuSaNWW0y9Fesb+ukgYG1qpM4KQRD5rSWcrh+hKOGPpyUl0I8Rlj9gpIZO/OMhWcoknlGi6cDzRDKI9ffAnEMQJjniwrgfberDvBBygoiFgEAQMBn/X3q97baFzOqCsAEpyoSQP7CUQEAKBoXofEAqBrn1AcP2/OmID63Io+R7s2yHodAiZCoEj3Z4rOEhiGYn8Sj/OpS9n4efBJZFCPj+ym3sZVWcXzk+G73u5DyuryoE8uOxdWPagzlXgLTYb3HKLXh03bVrT5Wkwk1dGWI8IRGDJEj3kj4mBBx/U8ab1YSmDWWvgszVw9jC4ZQL06VS7XaEDvi6Ab/N1nOW1UTCpXaVlKgjx5LKaNHaSSQyRnEpPjqNLhX80lwQOsYUMtgPQmRF040Q6MED7QJ1ZlSJq+wOcRVpAA0dD0GgIOAaUSyXKCiE7AbLjISte7/NTtJgq5co67Q/BkRAaBSFRWlDtpWAr1SJrK3E9LwFHmT43IAS6DIOux7q2YyDYjYkodrCuhNLvwbYFgsZD2C0QUHuRvYMydjOXfJIZw30e/b0PsYwnGV/x4+Pufzw+CX73YuLq8v/BPefAaYPqb1uOxaJXXF10kRbXpsS4Agytno0b4bHH9CTUzJnepYPLs8C7S+HnP+H6U2HJ/XqtfU2SrfBWDqwrgWsj4cfeeqlmOcXYWEYSK0mlL5GMozfXMwJ/V+WiYjKIJ5YDbKID/enGSAZyIYGE6Qs49kPpO1D6E6hwCBqrLcHw+8CvSnhBzj7Y+ymkrIGsPRDcHjoNgo6DoOdIOO4qiOoD/of5tbSVQuYuOLQddvwIX/0NHkyrPR5WARB8pt7EAZb3IPdy6Lyt1iX9CaYf57CC/5JPClHUVsYU8rHiIIdSetCu1vF8BzyaAWPqmd9KzYFp38OQbt6LqoieoHrhBf2DfOml3p3XVBhhNbQoiYlaUEXg3Xd1YuP6OJQPb/ymJ6NuPxN+f1ivTa/JWgu8maOD9O+Khle7Vh8W7iWXhewljULOpA9PM6HKML+MZNaQzHICCKYPZzKUqfiVf2Uch6D0UyidAyoUQi6HDnPBr8pEj9UC+xZC/K+Q/ocWzQHnwOkPazFtqjFqYAj0OFFvthLYv6l+J6MzE0o+h+j5bg8XkMIfvMUo7qklqoKwmH2sIIUHGEN3N6I6rxCez4KHO+lFAu4otcFrv8JvO+CFqXCql2Fza9fC44/DccfBggVNs3y1oRhhNbQImZnw7LM6k9Szz+q1+vWRU6SDxbelwb3nwvOX1dYLu8DsAngvFwYGwROdYXgVK9aKg99JZSlJdCaMSQxgMNGuWXwhmzj2sYQCUunFqYziXkJwJRFwFkDpd1D6I+APIZdBh+/Ar4M+LgIZu7SQJizWQ/R+E+H4a+CC1/WwvrlJXQe9x9bdRiyQdxW0fwf8a6fZzGQH25nFGB4gnK7VjlmwMYPNdCCEJ5lAYI36pJl2uOeg9qf+FqNzHNR6eYF5W+HFX+Dm8RD7oHc5XXfv1j/KoaE6EXZTZrJqKEZYDc2KxQJvvqkti0cf1Y/rM9wcTvhopZ4hfvRCeONqN3M8At8VaAt1cjv4thd0rvLpLsPOL+xlNWmcRi+mMZb2rlluwUEa60lgARH0YADnE8UAveJJBGzrwPIhOBIg5CqI+hz8XE5cEdi/RQ+5E2Oh42AYcgFc/gWEuTGdnE44uB+S90HKvsp9WkplkajyZKlRHaBLtypbVxhxAsT09f4NT1wG/c/wfFyckPcPCL8bgmq7C9NYzV5+ZRz/JYjqpmYiebzHZi5nKKfQo/plXf+PN3LghS46b4A7du2HB77TiWsW3qtji+vj4EF46im9KOTZZ+GEE+o/p7kxwmpoNpYu1WJ66606OYa/Fwbc6nh4+EeYdCysmOZ+dnhFMTyRqVdA1bSKynDwK3tZRSrn0p8XOYMAl1XloIwkYkkilm6cyBgeqGKd5kLJF1DyLQQeC+H/hsCT9DERSN2gxTRpJXQ7Do65DM58AvxrdHB/GvyxVm/bNutzu/eEPv2gd18YfRpcfi307F09el0E8nIh4yAcOgiZh7QIz3gD/nEbXHa1d296yhqY8JDn44WPQOAoCKntlIxnHpnsZByP4V8l1EoQFpLIWtJ4gDF0proaHrDBHQchJhCW9qm+rLecfIsefew6AK9eAcO9WKpaWKiTiS9bBv/9L5xzTv3ntBQmKsDQ5OTn6wkFiwXeeEMXy6uP/bkw7QetLy9O1csZa7KzDB7JgI7+8HRn6FlF06w4+I1ElpPC2fTlLPpVDFPLKGAvv3CAjfRhIn05mwBCXNbparC8D450CL0WQq4AP5fP8NBfsOkTLVY9RsKIy6DP+MrJJqcTtm6CDau1kO5PhR694OSxcMpYOPbEw68EWFoK0+6AdhHw5Mt1R71bLfD5RXDzUvfHLR/rCIbI96sNAQQn25iFEysncEu1CIBibPyPTXQhnGs4puJHCvTbNysfPsjV/uxT3VifTqcOh/tgBUybBJecVP+IpawMPvgAvvwS7rxTz/q3hkTfdUUFICJH9TZy5EgxtBxr1oiMHi3y00/etXc6Rd5bJjLueZGVu923sThE7j4gcmGyyLYSN68paXKfLJEFEi9WsVf83SFW2S6fyzJ5WFJklTiqHBPrFpGsM0Xy/k/Euq36BfdvFZl5jsi314jsXSbicFQ/brOJvP2KyLmjRR78t8icb0XSUry74cbgdIrcdIXI3bd4blNaKLLwAZGlT7k53yJS+LJI9gUiTmu1Q/mSKr/LM7JLfhCnOCv+7hCn/C6pcr8skY1yoNYl/ywRuSBZ5KGDIiWOWofF6RRZ/JfIhBdEnv5JxFJW/23a7SKffqo/P2+/LVLmxTnNCbBRPOhOiwtfU29GWFuOmTNFzj5b5NAh79oXWESumiEy7TsRq819mz9LRE5NFPkyr/axYrHKG7Je3pWNYpHqF8iVRImVabJXFlUTDHGWiRQ8IZI9ScS2r/oFbaUivz0m8ukkkexE9x3atkXkwgkiM97QStCUOJ0iSxbq13tqmkh2Vu021hKR1W+KzBgtsvkzEUeVPjmtIkUzRDJHixR/WE1Ui+SQbJA3ZZU8IzmSUHmKOOUP2S/TJFY+k21SKNXVbXepyNWpIlNTRHaUuu/2ijiRc14Vuf0zkfSc+m/T4RD5/nuRsWNFnntOpLCw/nNaAiOshmbFZhO56y6R227z3srYmS5y6nMiP29xf9zpFHkvR+TMfSL73Fxzt2TL/bJENkh69fPELrvkB1kuj0lRTUvLulkkc5xI0Xv6BaqSvFbkvVNF/vio9jERkZISkWcfEblykkiSB9H1JSlJItdeLPLA7SIHa1uM4nSKbPlCC+q6/4nYyqofs/wgkjlGpPB1EWelmV8iObJZ3pfl8phkyo5ql9wuGfKYLJcPZYvkSvWhQbJV5OZ0PWrYaHHf5TXxIue/JnLLJyJJmfXfotMpsnChyIQJIo88IpLjhQi3JEZYDc1GVpbI+eeLzJjh/TnfbRAZ/4LIXg+Wbb5dW0UPHhSx1tA4uzjlO9kpT8pKyZbq3/BC2S/L5TGJk9nirDrsd5aJFDzuslKTql+wrEhk/j0iX1wqkl9dpCtYu0rkvDEiX850L7q+xGYTefc1kcnjRTb/4b5N2iaRj8/S1nVZUfVjZb9rF0f+NBFHbuWfpVC2yacSK9PkgGyqZsXvkWx5SlbK2/KHZEhxtcsdtIncdUDkrCSRFTVeqpyN+0Qumi5yw0ciCV6MVkpLRb74QmTiRJG77xY5eLD+c1oDRlgNzcKOHSJjxoisWOFde6tN5N6vRf7xkWef25YSkTGJIvMLah/LEos8LivkR4kTRxVhcIpTEmShxMo0yZPkGi/qslKL368tiglLtcX35zfuBTM/T1uMN13u3mr0NVs36WH/O69qga1JUabInFtFvvybSPbe6sdscSI5l4nk3iBir3wPbGKRnfKtLJUHJFVWVxPUJMmT52S1vCrrJF2qj78zbCLTDoqM3yeysND927MlWeTSt0X+/r5I3P76by81VeSxx7QP9aWXRDK9sGpbE0ZYDU3OL7+InHaaSKKXo+L9uSJnvyLy/jL3X9KqQ/8Ua+3jayVN7pclskeyq/3dIlmySp6W7fKFOKTKiU6bSMHTItnn17ZSS/JF5vxL5Nu/ixRluO/Mj1+LnDNKZMEc727wcCgsFHn0XpG/XySSvK/2cbtNZM3bIu+NFdmzqMaxgyJ5t+uJKWulX8UhNomX+bJE7pN9sqSaBb9fCuVVWSfPy2rZJ9Wd14dsIg+4BHVuQe3/ldMpsjxOZPIbWlC3p9Z9a06nSGysyNSpIhdcIDJvXtO7ppsKI6yGJuXdd0UuvlgkP9+79uv3iox+RmSDBxEudYj8vY6h/3uyWd6UDbUmqA7IZlki90uWxNU4KV0ka6JI0XQRZ40p66TVIv8bJbJjrvvOFBdrP+oTDzT9LEpxscj3X+rogjnf1lYxh0Nk13yR908TWfVadT+qI0+k8DmRzNNESivF1ikOSZXfZYncL7tlrtir/NhkSLG8KxvlSVkpu2v8QB2yidx3QOT0fSLzPAjq/K0iZ7wk8n+feXbjlGO1inz4ocipp+rh/m4PER9HEnUJq1kgYDgsvvlGl7v48UfvAv53pMM9X8O8u3SBvpo4BG7Yr1dPXVejHLUgfMgWutGOS6ie9iqbOOL4gQk8VZkcBVzLNa+G9q/r1HhVsWTDwv/A9Qsg3ENwrcMBJSU6ZrQpKCyA3xbAgtmQmQHnTIbvfoXIKjdvydbxszt+hL4T4OofIMKV5dueDJbpYN0AYbdCxxWg/HFiJ4UVJLKILhzPeB6vWDmVQTE/EMdBipnKUI6rkj81yw6vZOuENQ92hFe61o4zXboTnpsPJ8TAV/9yX8amHKdTf0befhumToVFi6Bd7VQCRx1GWA2NZvVqnZl9/nzvRDX+ENz8CXx7m3tRFdErdsaE1hZVgC/ZQTiBtUS1kDS28jHjeKyGqDoh70YIu7u2qIrAT/8H577gWVRB10xuFwEHD0C37vXfpDfkZMOv82DBHCixaDF96lXo3ad6/1LXw/r/QUE6jLwRbl6mE6wA2DZB0Wsg+RB+D0S8BkrhwEoSi0gilp6MZTxPVrwnByniR3aT4RLUqgmp99tgeg6sscD9neDFLtUFVQRW7oHn58OALvD5P6FnB8+3KK5sUy+/rFdILVoE7esvRHDUYITV0CgSEuD++/WXx5uSwvsy4boP4Yt/us+VKgL3HYIeAXBPx9rHZ7ObYmzcSvWF4aXksoE3Gc19BFPDdCp6BgKOg9C/1b7gllk621Rd6+jLueASWPgT3Hhb/W09kZqsxXTxAv38/Ivhtfdri3VZEfz5Ffz5BXQeDuPug+7H62PihNIFUPwW+HeHdg9B4HEAWCkikUWks5Y+TGQiz1UsQ40nhznsxoqDSxjCCDpXvFxcGbyaDSk2XTn1hS6VJWhA1wn7YSPMWAbH9oL3rod+lae7ZdkyeOYZGDlSF/HzZqXdUYcnH0FLbMD5wDvdZl4AACAASURBVG4gAXjIzfF/AJnAVtd2S33XND5W35OVpWf/9+zxrn1qtsiYZ0V21TFT/OghvWrH3UTWQtkrr8v6ajP/IiJWKZZYeUiyxU1HLN+J5Fzj/oLZe0XeH6cXAHhDZobIlRd417Ycp1Nky0aRFx4XmTxO5MapIt9+5j6oX0QkI07k5zu1/3Ttu3pCreJaFh3FkDlWJP9BEXtaxaEiOSRb5EOJlWmSJLEVE3YOccoGSZfHZLm8KRskucak1JpikUtTRC5LEVlbPaJKRPRijdcXiYx6WuTZeSJZ9biXnU6RlStFJk0S+de/RNLS6m5/NMCR4GNVSvkD7wLnAGnAH0qpeSKys0bTb0XkjmbvoAHQS9Wvukpn+R/kRRLiA3lwxQz44AYY6mEk/VwmFDrhTTf+vBWksJWDPMCYanXrndhZz6sM5TKiqdER2yaw/A+if6l9QYcd5t4KF8/QpUy8oVNnsNsgP6+677MmWZmwdSMs+lknXDn2RDh/CtzzsHuz3umA3Qtgw3s6w//of8Pk6ZV9dhzS91G2CEKvgejfKvIW5JLAbuZgp4RBTOF4bkahsOFgKftYzD6G04l7GEVHV+VUp8CCIj3k7x+ordMhNd6CtBx4awmsTtAp/FY8BCF1pCMoLIQvvtDr+E84Ad56y/vy40czrUZYgVFAgogkAiilvgEuBmoKq6GFsNvh+uvhpptg3Lj622cWwNT/wdvX6GGkO97IhmQbvO+m3NMfHGA5yTzCqRXZ/EFPYm3iXXoyhu7UyIHhOAD5/4YOrgTUNVn1Mgy5ELqOqP8GqnLOZPjsQ7j9P9pHmpwIu/6CuL9g904oLYGOneGY4+HvN8FL73jOFFKcBZs+hp1zYOC5cMkHEFnlDbLtgOLXwZGss2q1exyULg1zkI3E8zMhdGAoUyuSThdh5VcSWUc6p9GLJxlPGFoRS53weT58kgenh8HnPaB7DbHcnAyvL4LsIrjrbF0Vta7kKDt36sTkf/4J110HCxdqd7RB05qEtSeQWuV5GjDaTbvLlFITgD3AvSKS6qaNoQl46ikYM0ZnF6oPEbh1ls5MNbKv+zbrLLC8GGb3rv0lLqCM79nFU0wgsEZ9pUNscZUKObf2RS3vQrv/ah+kO+J/hUmv138DNbnqBnj8frj0LG3B9u4LQ4+By6+DwcMgrJ5EorYSnfx62zdQnAEn3Qi3rKi0msUBZb+BZQaoMAj/DwTponwObKSynER+oxNDOYW7CEU7ogsoYx7xbCeD8xlQLS2iVbSYfpIHV7eH3/pAuxpavzMdHp+r0zE+eD4cXzvPdTXS0nTKvrw87WN/552mK4RwJNOahNXdv6dmTsOfga9FpEwpdRswCziz1oWUuhW4FSAmpp5PisErli2Dv/6C2bO9a//DRj3JMd5DMUC7wAOH4Ote4O/mP/8tO7mcYYS6+YjG8zMjud39ha0bIHya545d+S18eyVMeBiGTPLiTlxERsH0j7xvD7pY4O5fYOdsKDygrdMzH4fOQyvbOA5CyUwonQdBZ0L76RCgrdCqE1K9OI3xPFExw59LKXPZw26ymcIg/s4xFa4Su2gL9f1cuLI9LOsDoTUENTkLnvwJci3w9CX1l5c+cABeew02bdI/sBMmNOytaGu0JmFNA6r+e3sB+6s2EJHsKk8/BF5ydyER+QD4AHQ+Vt92s+2RlQUPP6yz/ntjneQUwWuLdHE/T7yVA5e1h15u/HdpFLKfIm6hW61jeewjiHaE4WZqWmwgJeBXx5i0fXcdt/rt1VrsTr6p/htqCJZsiPsZdswBa5GuJnDuCxDdv0o/nWBdpvO+SgGE3gwdV4LSia518cJ55JJAf87jDF6sqLWVTQlz2E0ieVzCYP7BsbrSAToG+JsCeCcHLolwn2Q6o0DHoO7cD09eXH+xvuRkHTIVFwf33ad968ZCrZ/WJKx/AIOUUv2AdOAq4O9VGyiluovIAdfTKcCu5u1i20NElxF++WXo6CYMyh0PfKe/tO08hGGl2mBuIcT2cX98Ftu4geMqBKMqu5nDYC5xf6JtS2WW/7oIjoBr5uhJrPxUbUU2Vi0cNkjbAPG/QfIq8A+GoRfpybH2PWq0zYCSz3TNrKDTIOJZCNAmvSDksJsE5mOlmMFcVDEhBZBEPj8TTwbFXMoQbub4imOFDpiZB1/lw+QIWOSmtlRyFry9FNYnwiOT4U035W2qsmePrnh68KBOUj5xohHUhtBqhFVE7EqpO4BFgD8wU0R2KKWeRoc1zAPuUkpNAexADjr8ytCEvPWWnu31dui3dCdYHXD+sZ7b3HtQZ5gPcPNF3cohIgmmb82YVKCQdGxY6ICHUq7WVRA03ruO+gfC32bC0idhzj+1ENYsq+IOEciO1/7ShMVQkgO9RsHAc2DCNAiq4WsVG5QtAMssoARCroWOy0DpXx0HZaSwkiRiaU8vhvC3igkpJ8JG9rOAvUQQxIUMZAiVv27JVng7B9aWwI1RsLxv7SH/2gR4czEUl8GdZ8MrV9QtkNu2aUEtKYGHHtI+dUPDMaVZDB7ZskVbKwsXQoAXP8GZBXDRW7DgHujoYdniz4XwaxG862ZuyYaTR1nOw5xKB2qbuxt4g4FMJhoPjtucSyHyPfDv6v64J/74UA/fr/gKgl0dt5VC7j7IToCcBNd+rx7eRw/QQjrg7MqlpbVuZjuUfALWtRA8CUJvgIBKE72IA+xlIdnE0Ztx9OFMglxloy3YiCWZFSRzHF24gIEVIVMA6y26SF+RE+6MhnPDq4ul3QGzN8F7y2FgF7jnnLprSonA4sV62WlYmHb7tMYCfa2NukqztBqL1dC6cDrhrrvgq6+8E1WA/3yrw3Q8iapN4NlMPTvtjliSGEdvt6JqpwwLmZ5FFcB5APzqWLjuiVP+CRE94NNzwc91s/7B2i/acSBED9QiGt0fgjyUGwVwFkHpbCj5FPx7Q+iNEPEqqEozMod44vgRwclAJnEc/0C5ZvGLsTGfeDZxkLPoy7OcTnCVr+g6CzyVqWt7/bcTHFPjbRLRgvrqIrjgWL10uHM9y0hXrtSz/KecooW1NZWQPpIxwmpwy7x5MH489K5ntricHel6uHn6EM9tZhfAxREQ6SGvwB8c4N+MdHvMQRlB1KMSodfpMtXhd3rX6aoMnaw3kYY5E52FUPYzlPwAkgvBU6DDD+BXWf3QiY101pHIb4QSzQj+Tnsqo1XyKGUe8fxFJhcwgBc4A/8q/uVNJVpQO/jDO91hQFC1HiDiSoyyAEb1g1/vhch6or+SkmDaNG2hfv019OhRd3tDwzDCaqiFCLz5Jnz/vffnPD0PHp9Sd5sPcuEbDwsFbDgowurWWgUtTn7Uk+kl7GbInqAnhryZxHKHN6LqLNDhUaU/6iQowRdC5JvgXz20r4QcElnEQTbRg1GM5r7K8troLFNz2EMK+VzEIK5lRLXVZdtK4clMCFU6y1TNVVIAaxLgqZ9gUFf48p/Qo47EKABFRfDii7B2rd6fckr9t2toOEZYDbVYvFj72DrXk2yjnG2pWoxPqCNkeFMJ9AmCzh4+cTvJ4hg8Z+tw4qgIOfKICoGoz6HwaZBiaPdIRZD9YSEC9p1gjYWyJSBFEHIRRL6lh/xVmyJks4sEFmCliP6cxzCuqPajkEoBPxJHPmVcwhBu5YRqERA7y7SFKuiy3iPc/NZsSdZxqNHh8P4N0LeeRCclJfDppzBrlnbxPPOMmeVvSoywGmrx6qswc6b37Z+aB09dXHebt3N09iRPbOIgo/E8HnVixw8vZu0DhkCHL8G+G4qeh6Knod3D2or1FhFwJLqENBYcSRBwDASf6Zocqz3zVkYBKawgjdVE0pdhXElkleG+IOwki5+JB+AyhjKI6m/IX6XwYjZYnPB4ZzjBjaDu2q8FVSntz/aUf6Gc/HyYMUNnmbruOoiNrX+RmOHwMcJqqMbKlXoCo5eHIXtNNiVBkD+MqKN9hh3S7XCim6X75cSRzQ14jtHSwupF0tdyAoZA1Cyw74WiF3QKwfCHIfj0yjZSBo5UcKTodfmOFHDsBXsC+PfXQhrxfMVKKHd9OsBGklmGnRJiOJ1xPF4tJ2wZdlaQQizJ9CWSaxhB7yq+YhFYWqwTo7Tzg/s7wkg379OWZHh5IZTa4ImL6x4dgI4/nT5dJyG/7Ta9D/Tid8ngG4ywGioQ0TPEX37p/TkvLIBnLq27zcw8uKWOpFA5lBBFSLVEKzXxwx8bJd53rJyAARD1EdiToPglKHwUlD96oB2kh/L+fbR/NGgc+F8H/v3qHCdbyCSe+WSxg26M5HhuJLzGKrFCrMxlN9vI4HRieJxxFUlRQL/XC4rg5WwYGQJvd4O+QTVfSbtZHvkROoTDIxd6TmZTTn4+PPEEbN8O994Lzz3nOReMoekwwmqoYNEiGDXKe2s1zwL5JTCsnhnlVRa4rY44SitOwur5KLajOxYyKCWXEOqZoXFHQF+InAHOYp31SjVMbZzYOcRWklmGjRIGcSHHcUNFqFQ5xdiYxx62cIhLGMw1NSakALaUwMMZMDQY5vSCjm5uPadIJ0dJyYHXr4LBHsJlK/rn1On7/vc/HYf6xhvGh9qSGGE1VDB9Onz4offt5/8JFx5fdxsRyHNAVB2j+GD8KcNR53UUfgznSnbyDSfxf953siZ+dcSh1kBwksUuUlhBHvvoyvEM5yraUzsGLZsSfiGB7WQymYFcwfBqIVMAiVZ4OlMH9r/dDQa5meV3OOGDFfDZGnjsQphcz/sLepj/3//CGWfA0qUQ7v0tGpoII6wGADZv1lEA3lqroIPRp/+97jb7bNDPzRC3KsH4Y61HWAG6ciIJLGQ3cxjIZPyp58KNQBDySCSVlWSyg44MpR9n04FBbnMXJJLHT+whj1ImM9CthRpXBi9kQZYDHukEp7mZPBKBX7fD8wvgohNg+YM6lV9dJCZq6zQ0VLtvetYxKjA0L4ctrEqpaSLiNsuU4cjhtdd0wLi3FJVqV0DvOmb6Af4ogVH11MTyxmItZyzT2MdvLOcRBjGF3ox3K3gNwUYxGWwng23kEE8kfYhhPMdyPcrNhJkTYRMHmU88kQRzMYMZ4MY9sb0Uns+CUoGHO8EoN5NSIrB4B7z4CxzfG76/ve6qp6D9qM8/r38MX3gBTna7qNLQkjRYWJVS31V9CpyAh/R9hiOD5GSduPi447w/Z+F2mFRHopVy/ijR6QHrwh8/nLVS77rHD38GMIkYTmc3s1nBYxzD3+nMMV6dDzpxdDZxZPAnWezCjwA6M4I+nMEJ3OxWTEGv4V9BCstI5hg6cxenVFvDX86mEi2o/kpbqO7CpkQgdpee/BveQxdZrC+4326Hjz+GTz7RSaZffNH4UVsrjbFYC0TklvInSqkZPuyPoQWYPl3PIDeE2ZvgOTfFT2uytQye8aKKqx1ng14/kDBGcC0WsviLL4jnZwZwHh0ZihM7JWRjIZtSsikljzLyKaOAUnIRnHRkCF04gWFcXlHN1BNxZLOQvWRiYRy9eYoJbhNwr7XoIX97f3imCwz3cNl1e/XE1MAu8OnN0Kseq19EJxh/9VW49FJYvty7yriGlqMxwvpcjeeP+qIjhpbB4dDLG197zftzRCA1B/p3qb+tTWqnsnNHNCEcoIjueMjg4oEwOjGKe8hjH+msZTdz8SeIUDoSSjShdKQDAwgmimAiCSYSfy8WGgjCZg7xE3voQhh/Ywh93KQyBIgvg0cy9H2+3g0GenD95hbDY7MhoxA+vMF9GfCa7N8Pd94JMTHw668Q2YgcM4bmp15hVUr1Bf4NDEDnQN2qlPpZRJIBRCSnKTtoaFrWr9c5NxsypNyfBz29jHjy9rJTGMxsdntMwlIfUfSryGN6ODgQ1pDGAhIYRDR3exjuA2Ta9Sx/vBWe7wIneVgAYXfAhyv1TP8jk/XkVH2I6CH/xx9rS3Xs2MO4KUOz400w309AHJWlqY8HViql3lVKeVk/2NBa+eknuLie5ag12ZgEJ/f1bT+G04mDFJHdmEUAPiATC9+xi2nEkkYBj3AqN3O8W1EtccKLWXBJKkxqBwtjPIvqsl0w8SU92bf8Qe9ENSkJLrxQz/ovXWpE9UjEG1eAv4h8DKCUyhGRfyqlAoB70XWlbmjKDhqallWr9OqchrBxH5wxzPd9uYTBzGMPN+JF8KYPKMPBOtJZTjL++HEmfbiUibWqwpZjF/gyH2bk6pVkK/q6r4IAsC8THvwe2ofAD/+uf6YftFtmxgz47jvt9z7xxEbfmqGF8UZYlyil7hCRd3BVTRURO/CKUmpPk/bO0KTExcGgQd4nsi5nUzLcd77v+3MS3fieODKx0JmmyRQiCHvIYSlJJJHPWHpyJycT7WG4D9pC/TQPPsuHi9q5L9JXTlGpnulflwgvTYWTvfROrFoFjz4KF1ygrVSzrv/Ixpuv1H+Ah5VSG4EertLSFmAskF3nmYZWzdy5DXcDiOhlrFFNoHsKxf9xEi+zlhs5nuF1pBFsCPmUsZ0MtpNJIrn0pwPn0I+BdKgzBjbfoa3TuYVwXaQW1DAPgmqzw8zfYeYquPscePZv3vmtU1N1/HBAgE44bYL8jw7qFVYRcQLPKaXeAM5Gx612AP7CRAQc0cTGwh13NOyc7CLoUkd16ZoIYBUI8nIWqw+RPMY4PmIrX7OD8cQwjl7VEph41U9KWEsa69lPIH4cT1fOoz99iKy11LQmToFP8uCjPLi9A6zqC4F1nLI5Ge75GqacAMunQagXC8JE4PPP4YMPdElp40c9uvB6ECgiFmCeazMc4djtUFoK7RoW3cTBfOhRR6aqmpwSohcJuFvG6YlIgrmP0RRiZSUpPMtqOhPGWfRlBJ0rlozacJLvilLNo5Q8SsmihO1kEEYgY+nJNMbSrgFLX9dZ4KEMOCcclvWBkDqmd0usOjfqrgM6wD/Gy/Lgublw++3QvTssWWJiUo9GTK6ANsqOHTBiRMPPO5jv3URMOWeF63yjDRHWciIIYjIDmcxAksgnliQ+YzuBrmCWAPyIJJgoQlxbMIOIZgqDGmzhrrfAc1k6WcznPaF3Paev2A0P/QD/PhNenOp9uNqyZXp9/zPPwDnnNKiLhiMII6xtlHXrGlcz/mBBw4R1fDi87oNI575EclMTRAussejlp5384dWuMLieAMKCEpj2vfYz/3QndKlnuW45RUXwyCNw6BAsWAAdvbRuDUcmJgVuG2Xdusb59fbnQfcGCGs7P+2zLG7YitUmZ1UxXJACH+fB9G7wac+6RdXphC/Wwrmv6RwJX/3Le1FdtAjOPhtOPx2+/daIalvAWKxtlPh4GDiw4eel5cB53uc7AWBcGPxugfMa6M/1NSIQW6yz9vcJhHe71Z/SEHQl1Ednw4TBsPQBCPdyWUxODvznP9pN8MsvEF1PTgDD0UOrElal1PnAdMAf+EhEXqxxPBj4DBiJDvW6UkSSmrufRzoFBdC+feMyIx3Mh64NXK8+JULXdGopYbW5Avs/yIXRofBBd10xtj4KSuDer6HMDp/dUn+KxKqU+1KffhrOPbfxfTccmbQaYVVK+VO5bDYN+EMpNU9EdlZpdjOQKyIDlVJXodMVXtn8vT2yiYuD4cMbd65S3q//L2dkKJSJtlrHNWOFUKvA566wqUsjYFEMRHhZj3Dlbr1y6rEL4UIvlqGW43TqdH7r1xsrtS3Tmnyso4AEEUkUESvwDVAzfP1iYJbr8Q/AWUqZjJQNZdcuGNbIJakB/mBvhL/0za4w7ZAWu6amzAnv58KEJMhzwpI+8GAn70S1zKYnp6YvgZ/vapioZmXB31wLA+bMMaLalmlNwtoTSK3yPM31N7dtXMtq8wEzFdBA4uJg6NDGnRvgBzbvkv1Xo3sgXB8FL2c17nW9IdMO07Ph9GS9DDW2D9zX0fPy05qsTYCzXoGh3eCH26Gzl5NToIf+kyfDPfdoF4CpjNq2aTWuANyPMGvaN960wbXs9laAmJh6CrC3QeLidAb6xhDor9PgNYZ/RsG5KTCsQA/N/Xww1ihywtwC+KZAJ0m5or0O7PcmB2w5WYU6JrWoDL79P+9TIgLYbPDkk/o9/eUXM+Nv0LSm39U0qFb+shew31MbV4atSHSO2GqIyAcicrKInNy5c+cm6u6RS3Z24wWgsa4A0EL6XS9duuT0JJhdoEOxGopVYH4hXJ0Gk1Mg06EnpH7tAzd18F5UnU74YDlc9BZccQp8c1vDRHXfPjjvPL2+/4cfjKgaKmlNFusfwCClVD8gHbgKqFkDdB46TeFaYCoQKyLN4LU7erDZGp7NqipB/toP2Vii/eH5rpBt1wsH3siGKyN11ih3M/U2gb1W2FXm2qyQYIXzwuHJzjCkkRmB9xyE2z6DM4d5VxG1JrNnw+uvw3vvNW4Fm+HoptUIq4jYlVJ3AIvQ4VYzRWSHUuppYKOIzAM+Bj5XSiWgLdWrWq7HRyYHDkCPHo0/v0t7vfrqcNdAdQyA57pAnkNbnw9kwAEbnB6urdg4q/aZBioYEATDguHEELgqEvoFHl4Rvc/X6Iz+H/0DBndr2LkiOn/tX3/Bb79BWDNGORiOHFqNsAKIyC/ALzX+9niVx6XA5c3dr6OJ9PTDS003qKu29s7zkZUW5Q/XRumtzKlDsoIU3BsMnX386cwthnu/gYgQ+O0+CGmglWqxwD//CUOG6BR/Jh7F4IlWJayGpicpSRemaywjesLbcT7rTjWC/eCsJlhEIALf/QFv/AZPTIFJDSjzXU5SElx/Pdx9N1x2mc+7aDjKMMLaxti7F045pfHnl1usRwop2XDnl7rUdEOWo1YlNlZn9//wQ+NPNXiHEdY2xt69cNVheKaDAnQcq0jrHgo7nPDOUvhxE7xxFYzs2/BriMCbb2ph/eUX6NCAiAFD26Y1hVsZmoGUFOjT5/Cu0bMDpOf6pj9NwV9pOtDfaofYBxonqqWlcMMNkJmpS9gYUTU0BGOxtjHs9sMvVHfGUJj/J9x2hm/65CtEYNZqmLUGZt4I/bs07jp5eXDFFXDTTYdn3RvaLsZiNTSYq0fDV+tbuhfVScqCKW/pMim/3tt4UU1Ph4sugoceMqJqaDzGYm1DHO7igHLah8LgrrApqXHDbF9id8BbS+CnLfDm1XDiYbg5duyAW26BGTPghAYkXzEYamIs1jZETo7vMi7dejp8sMI312osW1Pg7Ff146UPHJ6oLlsG//oXfPONEVXD4WMs1jZEdjZ06uSba53SD+77FopKoV0zVxktscJTP8H2dPjkJuh3mOkgPv8cvvwS5s+HqAZUoDUYPGEs1jbE4SRfqYlScNUo+GaDb67nDSKwcBuc+QqM6AXz7z48URXR1VIXL4Z584yoGnyHEdY2RG6ub8XjmrHw0UrIKfLdNT2xNUUX8vv1L5h3J1w79vDiaJ1OuOsuHVY1axYEeVGqxWDwFuMKaEOUlUGID4ftUWHw0uVw/Ufw013g3wQ/01mF8NgcOJQPM66DgV0P/5pOJ9x2G/Trp5NSGwy+xlisbYiyMghuZJo9T5w+BM4eDv+d49vr2h3w9hIdQnXpSTDnTt+IqsOh41OHDjWiamg6jLC2IZpCWAHuPkevyZ+9yTfXW7oTJr6kE2qvmOa7TFo2m06kcsopuiy1wdBUGFdAG6KsDNo1QfYopeCDG+C812FYdxjWiHyvh/Jh3la9tr9XB5h9h8796issFrj6arjkErjxRt9d12BwhxHWNkRTWawAYcEw62a49kP4+EadQLo+n2tSFszZDAv+hNAgmHKCvkbXSN/2LT8fpk6F22+HSy/17bUNBncYYW1DWK1NO/vdvwu8eiW8tgjiD+kqjz2iYEBnGNBFb307wrI4+HQ1dInQ/tMf/w2RTZSJPzdXW6mPPw5nndU0r2Ew1MQIaxvCz0/HbjYlpw7UG+jZ9wP5sDdDb4t3QGImnNwXfrqz6cS0nIIC+Nvf4Omn4fTTm/a1DIaqGGFtQ/j76+xWzYWfn04x2LMDTBjSfK8L2qc6dapOpmJE1dDcmKiANkRAgA43OtopLYXLL4c77tDlqQ2G5sYIaxsiIKB5LdaWwGrV6f6uvx6mTGnp3hjaKkZY2xDN7QpobqxWuOYaPfN/5ZUt3RtDW8YIaxviaLZYrVYtphdeqEuqGAwtiRHWNkRYmJ7UOdpwOnWC6osvNqJqaB20CmFVSkUrpRYrpeJde7el25RSDqXUVtc2r7n7eaQTEQGFhS3dC9/z8MMwbBj84x8t3RODQdMqhBV4CFgqIoOApa7n7igRkRNcm5maaCDt2+vYzqOJN9/UVvhDnj4xBkML0FqE9WJgluvxLOCSFuzLUcvRZrF+8w2sXavF9XBysxoMvqa1CGtXETkA4Np7qrEZopTaqJRap5Qy4ttAIiKOHos1NhZmztRJqv39W7o3BkN1mm3llVJqCdDNzaFHG3CZGBHZr5TqD8QqpbaLyF43r3UrcCtATExMo/p7NNK+/f+3d/+xVpR3HsffHyHYrKURAqsoiGilu+BW0KO1Wk1s0WI3rRVps9a1uiCWXbfbrtvUNm5Ms7Zdsyau2pqwuH8ojdW1GKJAlYKRtdYYuShUqEv5ZRQwCq6tXZQay3f/mKEc8dzrufc+M3POnc8rOblz5sw53+fRuR/mPnPmmaFxxLp+fXbt/7JlaSfuNkultGCNiBm9vSbpFUnjIuJlSeOAV3v5jF35z22SVgPTgfcEa0QsBBYCNBqNgq+O7x5HHplNStLNXn01u5vq4sW+R5V1rk4ZCngIOPBFmSuABw/dQNIoSYfny2OAs4FfldbCIWDEiGyy5261b192AcAtt8D48VW3xqx3nRKsNwHnS9oMnJ8/R1JD0n/m2/w50CNpPfAYcFNEOFj7qejZrYqyf382QfXcuXDWWVW3xqxvHTG7VUS8BrxntsyI6AGuypefBP6i5KYNOcOGZROxdNsJn29+E049NZsHJB3BuQAAC3pJREFUwKzTdcoRq5VkzBjYs6fqVvTPD36QDQN84xtVt8SsPR1xxGrlOeooeOWV7Gc3WLIEVq+G++/3d1Wte/iItWaOPhpefrnqVrRn3Tq47Tb40Y+6b+jC6s1HrDUzYQLs2FF1K97fb34D8+fDT36STR5j1k18xFozEyfC9u1Vt6JvB74B8N3vZv8QmHUbB2vNTJ4MmzdX3Yq+3XQTnH46zOj1khKzzuahgJoZN66zx1hXrYI1a+CBB6puidnAOVhr5sCZ9YjOO8ve0wPf+Q4sXZrd4dWsW3n3raFjjoGdO6tuxbtt2QJf/Wp2smpUy2nOzbqHg7WGpk+HZ5+tuhUH7d2bzf5/113ZUIVZt3Ow1lCjkY1jdoIImDcvu6rqIx+pujVmaThYa+i002Dt2qpbkbn1VjjhBPi8py23IcQnr2royCPht7+t/gTW6tXw6KPw4HsmiTTrbj5iramJE+HFF6ur/9JL2Q0AFy3y5ao29DhYa6rRgKefrqb222/Dl78MCxbA6NHVtMGsSA7Wmjr33OxP8SrceCNceilMm1ZNfbOiOVhravp0eOaZ8u8osGZNdjPAefPKrWtWJgdrTR12GHz4w9kX88uydWt2EcCCBZ131ZdZSg7WGps5Ex55pJxaO3dmNwJctCi78stsKHOw1tgFF8DKlcXX2bMHvvjF7Eh18uTi65lVzcFaY2PHZt9n3bevuBpvvAGzZ8PNN/tkldWHg7XmPvUpWLGimM+OgDlz4LrrfMtqqxcHa81dfnk27lmE22+HqVPhwguL+XyzTuVLWmtu0iT43e9g9+5saCCVp5+Ghx+G5cvTfaZZt/ARq3HZZXDvvek+7/XX4Wtfy6YB9OWqVkcdEaySviBpo6T9khp9bDdT0iZJWyR9q8w2DmWXXAKLF6f5rAiYOxe+973sVttmddQRwQpsAGYBj/e2gaRhwB3AhcAU4FJJU8pp3tD2wQ/CiSdmV0QNxltvwde/DqecAp/8ZJq2mXWjjgjWiHg+Ija9z2ZnAFsiYltEvA3cB1xUfOvqYc4cuPPOgb//xRfh/POzyV1uuCFdu8y6UUcEa5uOBV5qer4jX2cJfOITsGFD9mX+/lq37uAFAJdf7stVzUoLVkmrJG1o8Wj3qLPVr2vLKUQkXS2pR1LP7t27B97oGpHgmmvgjjv6976VK7Pr/x94AE4+uZi2mXWb0oI1ImZExMktHu3OH78DmND0fDywq5daCyOiERGNsSm/QzTEzZqVzR2wd2972y9aBLfcAsuWwbH+28Hsj7rpe6xrgJMkTQJ2An8FfKnaJg0tw4bBV76STc4yYwacfTZ87GMwcuS7t4uA738fNm3KbqsyYkQ17TXrVB0xxirpYkk7gI8DyyWtyNcfI+mnABHxDvD3wArgeeD+iNhYVZuHqiuvhKVL4cwz4ec/z67zP+ecbJjgnnuyqf/mz4c334S773aomrWiKHum45I1Go3o6empuhld7Z13shNbv/gFPPkknHceXHVV1a0yq5aktRHR8nv33TQUYBUZPjybmWratOzI1cz61hFDAWZmQ4mD1cwsMQermVliDlYzs8QcrGZmiTlYzcwSc7CamSXmYDUzS8zBamaWmIPVzCwxB6uZWWIOVjOzxBysZmaJOVjNzBJzsJqZJeZgNTNLzMFqZpaYg9XMLDEHq5lZYg5WM7PEHKxmZok5WM3MEnOwmpkl5mA1M0usI4JV0hckbZS0X1Kjj+1ekPScpHWSespso5lZu4ZX3YDcBmAW8B9tbHteROwpuD1mZgPWEcEaEc8DSKq6KWZmg9YRQwH9EMDPJK2VdHXVjTEza6W0I1ZJq4CjW7x0fUQ82ObHnB0RuyT9KbBS0v9ExOMtal0NXA1w3HHHDbjNZmYDUVqwRsSMBJ+xK//5qqQlwBnAe4I1IhYCCwEajUYMtq6ZWX90zVCApCMkjTywDFxAdtLLzKyjdESwSrpY0g7g48BySSvy9cdI+mm+2VHAE5LWA08DyyPikWpabGbWu075VsASYEmL9buAz+TL24BTSm6amVm/dcQRq5nZUOJgNTNLzMFqZpaYg9XMLDEHq5lZYg5WM7PEHKxmZok5WM3MEnOwmpkl5mA1M0vMwWpmlpiD1cwsMQermVliDlYzs8QcrGZmiSliaN+5RNJuYC9Q5S2zx1RY37Vduy71y649MSLGtnphyAcrgKSeiGjUsb5ru3Zd6lfd92YeCjAzS8zBamaWWF2CdWGN67u2a9elftV9/6NajLGamZWpLkesZmal6fpglTRT0iZJWyR9q8Xr50p6RtI7kmYf8tofJK3LHw8VUPtaSb+S9EtJj0qa2PTaFZI2548rSq49qH63WX++pOfyGk9ImtL02rfz922S9Omyaks6XtJbTX1fkLp203azJYWkRtO6QvvdW+0y+i3pSkm7m2pc1fRa0ft6X7UHva8PSER07QMYBmwFTgBGAOuBKYdsczzwUWARMPuQ1/6v4NrnAX+SL/8t8F/58mhgW/5zVL48qozag+13P+p/qGn5c8Aj+fKUfPvDgUn55wwrqfbxwIYi+51vNxJ4HHgKaJTV7z5qF95v4Erghy3eW8a+3rJ2in19oI9uP2I9A9gSEdsi4m3gPuCi5g0i4oWI+CWwv4Laj0XEm/nTp4Dx+fKngZUR8b8R8TqwEphZUu0U2qn/RtPTI4ADg/kXAfdFxO8jYjuwJf+8MmoP1vvWzt0I/Buwr2ld4f3uo/ZgtVu7lcL39U7U7cF6LPBS0/Md+bp2fUBSj6SnJH2+4NpzgYcH+N6UtWFw/W67vqRrJG0l+0X/hwG2PWVtgEmSnpX035LO6UfdtmpLmg5MiIhlA2l3QbWh4H7nLsmHnhZLmtDP9xZRGwa/rw9ItwerWqzrz9HJcZFdqfEl4FZJJxZRW9JfAw3g5v6+t4DaMLh+t10/Iu6IiBOB64B/7s97C6r9MlnfpwPXAj+W9KFUtSUdBvw78E8DbXdBtQvtd24pcHxEfBRYBdzdj/cWVRsGv68PSLcH6w6g+V+n8cCudt8cEbvyn9uA1cD01LUlzQCuBz4XEb9P0e5B1h5sv9uu3+Q+4MDRQil9b1U7/zP8tXx5LdnY3eSEtUcCJwOrJb0AnAk8lJ9EKrrfvdYuod9ExGtN+9idwGntvrfA2in29YGpYmA31QMYTjYYPomDA9tTe9n2LppOXpENpB+eL48BNtPiZMBgauf/E7cCJx2yfjSwPW/DqHx5dEm1B9XvftQ/qWn5s0BPvjyVd5/E2Ub/TuIMpvbYA7XITobsTP3f/ZDtV3PwBFLh/e6jduH9BsY1LV8MPFXivt5b7UHv6wN9FF6g8A7AZ4Bf5yFyfb7uX8iO0gBOJ/tXby/wGrAxX38W8Fz+P+o5YG4BtVcBrwDr8sdDTe+dQ3YCYwvwN2XVTtHvNuvfBmzMaz/W/MtAdhS9FdgEXFhWbeCSfP164Bngs6lrH7LtavJwK6PfvdUuo9/AvzbVeAz4sxL39Za1U+3rA3n4yiszs8S6fYzVzKzjOFjNzBJzsJqZJeZgNTNLzMFqZpaYg9XMLDEHq5lZYg5WqzVJwyTdJmljPofrCVW3ybqfg9Xq7tvAtoiYCtwO/F3F7bEhYHjVDTCriqQjgIsj4sCkHduBv6ywSTZEOFitzmYAEySty5+PJptjwWxQPBRgdTYNuCEipkXENOBnwDpJR0i6W9Kdki6ruI3WhRysVmejgDcBJA0HLiCbNHkWsDgi5pHdM8usXxysVme/JpsQGuAfgeWR3Y9qPAdvB/KHKhpm3c3BanV2L3CqpC1kd/K9Nl+/g4M3X/TviPWb52M1O0T+bYEfkt3p9ImIuKfiJlmXcbCamSXmP3PMzBJzsJqZJeZgNTNLzMFqZpaYg9XMLDEHq5lZYg5WM7PEHKxmZok5WM3MEvt/BKrSCCePAR4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "DIR = 'results/OP'                                             \n",
    "RERUN = not utils_os.is_file_exist(DIR, 'true_samples.npy') \n",
    "\n",
    "## Define the problem\n",
    "problem = problem_OP.OP_Problem(N=2000, n=50)\n",
    "true_theta = problem.get_true_theta()\n",
    "\n",
    "## Get x_o ~ p(x|theta)\n",
    "if RERUN:\n",
    "    # observed data x_o\n",
    "    problem.data_obs = problem.simulator(true_theta)\n",
    "    problem.y_obs = problem.statistics(data=problem.data_obs, theta=true_theta)\n",
    "    utils_os.save_object(DIR, 'data_obs', problem.data_obs)\n",
    "    utils_os.save_object(DIR, 'y_obs', problem.y_obs)\n",
    "    \n",
    "    # true samples theta ~ pi(theta|x_o)\n",
    "    true_samples = problem.sample_from_true_posterior()\n",
    "    utils_os.save_object(DIR, 'true_samples', true_samples)\n",
    "else:\n",
    "    # load previously simulated true samples & x_o\n",
    "    true_samples = utils_os.load_object(DIR, 'true_samples.npy')\n",
    "    problem.data_obs  = utils_os.load_object(DIR, 'data_obs.npy')\n",
    "    problem.y_obs  = utils_os.load_object(DIR, 'y_obs.npy')\n",
    "    \n",
    "## Visualize\n",
    "problem.visualize()   \n",
    "visualization.plot_likelihood(samples=true_samples, log_likelihood_function=problem.log_likelihood, dimensions=(0,1))\n",
    "plt.savefig('OP_true_posterior.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SMC-ABC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Sequential Monte Carlo ABC\n",
    "\n",
    "hyperparams = ABC_algorithms.Hyperparams()\n",
    "hyperparams.save_dir = DIR\n",
    "hyperparams.device = 'cuda:1'\n",
    "hyperparams.num_sim = 4000                        # number of simulations\n",
    "hyperparams.num_samples = 200                     # number of samples to represent posterior\n",
    "hyperparams.L = 2                                 # number of rounds in sequential learning\n",
    "\n",
    "smc_abc = SMCABC.SMC_ABC(problem, discrepancy=discrepancy.eculidean_dist, hyperparams=hyperparams)\n",
    "smc_abc.run()\n",
    "\n",
    "JSD_smc_array = []\n",
    "for l in range(hyperparams.L):\n",
    "    print('round =', l)\n",
    "    smc_abc.posterior = smc_abc.posterior_array[l]\n",
    "    visualization.plot_likelihood(samples=true_samples, log_likelihood_function=smc_abc.log_likelihood, dimensions=(0,1))\n",
    "    JSD = discrepancy.JSD(problem.log_likelihood, smc_abc.log_likelihood, true_samples, true_samples, N_grid=30)\n",
    "    JSD_smc_array.append(JSD)\n",
    "    print('JSD smc = ', JSD)\n",
    "utils_os.save_object(DIR, 'JSD_SMC', JSD_smc_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Sequential Monte Carlo ABC +\n",
    "\n",
    "hyperparams = ABC_algorithms.Hyperparams()\n",
    "hyperparams.save_dir = DIR\n",
    "hyperparams.device = 'cuda:1'\n",
    "hyperparams.num_sim = 4000                       # number of simulations\n",
    "hyperparams.num_samples = 200                    # number of samples to represent posterior\n",
    "hyperparams.L = 2                                # number of learning rounds\n",
    "hyperparams.type = 'cnn1d'                       # the network architecture of S(x)\n",
    "hyperparams.stat = 'infomax'                     # statistics function: infomax/moment/score   \n",
    "hyperparams.estimator = 'JSD'                    # MI estimator; JSD (accurate) or DC (fast)\n",
    "\n",
    "smc2_abc = SMC2ABC.SMC2_ABC(problem, discrepancy=discrepancy.eculidean_dist, hyperparams=hyperparams)\n",
    "smc2_abc.run()\n",
    "\n",
    "JSD_smc2_array = []\n",
    "for l in range(len(smc2_abc.posterior_array)):\n",
    "    print('l=', l)\n",
    "    smc2_abc.l = l\n",
    "    smc2_abc.posterior = smc2_abc.posterior_array[l]\n",
    "    visualization.plot_likelihood(samples=true_samples, log_likelihood_function=smc2_abc.log_likelihood, dimensions=(0,1))\n",
    "    JSD = discrepancy.JSD(problem.log_likelihood, smc2_abc.log_likelihood, true_samples, true_samples, N_grid=30)\n",
    "    JSD_smc2_array.append(JSD)\n",
    "    print('JSD smc2 = ', JSD)\n",
    "utils_os.save_object(DIR, 'JSD_SMC2', JSD_smc2_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SNL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Sequential Neural Likelihood\n",
    "hyperparams = ABC_algorithms.Hyperparams()\n",
    "hyperparams.save_dir = DIR\n",
    "hyperparams.device = 'cuda:1'\n",
    "hyperparams.num_sim = 4000\n",
    "hyperparams.L = 2\n",
    "\n",
    "print('\\n SNL ABC')\n",
    "snl_abc = SNLABC.SNL_ABC(problem, discrepancy=discrepancy.eculidean_dist, hyperparams=hyperparams)\n",
    "snl_abc.run()\n",
    "\n",
    "JSD_array = []\n",
    "for l in range(len(snl_abc.nde_array)):\n",
    "    print('l=', l)\n",
    "    snl_abc.nde_net = snl_abc.nde_array[l]\n",
    "    visualization.plot_likelihood(samples=true_samples, log_likelihood_function=snl_abc.log_likelihood, dimensions=(0,1))\n",
    "    JSD = discrepancy.JSD(problem.log_likelihood, snl_abc.log_likelihood, true_samples, true_samples, N_grid=30)\n",
    "    JSD_array.append(JSD)\n",
    "    print('JSD snl = ', JSD)\n",
    "utils_os.save_object(DIR, 'JSD_SNL', JSD_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Sequential Neural Likelihood + \n",
    "hyperparams = ABC_algorithms.Hyperparams()\n",
    "hyperparams.save_dir = DIR\n",
    "hyperparams.device = 'cuda:1'\n",
    "hyperparams.num_sim = 4000                       # number of simulations\n",
    "hyperparams.L = 2                                # number of learning rounds\n",
    "hyperparams.type = 'cnn1d'                       # the network architecture of S(x)\n",
    "hyperparams.stat = 'infomax'                     # statistics function: infomax/moment/score   \n",
    "hyperparams.nde = 'MAF'                          # nde; MAF (D>1) or MDN (D=1)\n",
    "\n",
    "snl2_abc = SNL2ABC.SNL2_ABC(problem, discrepancy=discrepancy.eculidean_dist, hyperparams=hyperparams)\n",
    "snl2_abc.run()\n",
    "\n",
    "JSD_array = []\n",
    "for l in range(len(snl2_abc.nde_array)):\n",
    "    print('l=', l)\n",
    "    snl2_abc.set(l=l)\n",
    "    visualization.plot_likelihood(samples=true_samples, log_likelihood_function=snl2_abc.log_likelihood, dimensions=(0,1))\n",
    "    JSD = discrepancy.JSD(problem.log_likelihood, snl2_abc.log_likelihood, true_samples, true_samples, N_grid=30)\n",
    "    JSD_array.append(JSD)\n",
    "    print('JSD snl+ = ', JSD)\n",
    "utils_os.save_object(DIR, 'JSD_SNL2', JSD_array)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
